#include <stdlib.h>
#include <cmath>
#include <vector>
#include <iostream>
using namespace std;

//分段函数多项式类
class PiecewisePoly
{
 public:
  vector<double> coe;//系数
  vector<int> exp;//次方
  double begin;//区间左端点
  double end;//区间右端点

  /* 
  void dis(){
    for(int i=0;i<coe.size();i++){
      cout<<coe[i]<<"x^"<<exp[i]<<"+";
    }
    cout<<endl;
  }
  */
  
};

class Solve_Base//求解基类
{
 public:
  vector<PiecewisePoly> answer;
  vector<PiecewisePoly> answer_;
  vector<double> t;
  int N;
  vector<double> x;
  vector<double> y;
  
  Solve_Base(int _N)
    {
      N = _N;
      // input();
    }
  Solve_Base(int _N,vector<double> _x,vector<double> _y)
  {
    N=_N;
    if (_x.size()!=_N || _y.size()!=_N)
      {
	 
     cout << "Error Input!!" << endl;

     exit(0);

      }
    for (int i=0;i<=N-1;i++)
      {
	x.push_back(_x[i]);
	y.push_back(_y[i]);
      }
  }
  /*
  
  void input()
  {
    int i;
    double xin;
    double yin;
    cout<<"input "<<N<<"numbers"<<endl;
    for (i=1;i<=N;i++)
      {
	cout<<"x"<<i<<"= ";
	cin >> xin ;
	cout<<"y"<<i<<"= ";
        cin >> yin ;
	x.push_back(xin);
	y.push_back(yin);
      }
    
  }
  */
  


};


class ppForm : public Solve_Base
{
 public:
  ppForm(int _N)
    :Solve_Base(_N)
  {
    
  }

  ppForm(int _N,vector<double> _x,vector<double> _y)
    :Solve_Base(_N,_x,_y)
  {
    N=_N;
    for (int i=0;i<=N-1;i++)
      {
	x.push_back(_x[i]);
	y.push_back(_y[i]);
      }
  }
  vector<PiecewisePoly> S_01();
  vector<PiecewisePoly> S_23(int,double,double);
  int condition;
  double boundA;
  double boundB;
  
};

vector<PiecewisePoly> ppForm::S_01()
{
  int i;
  double coe_0;
  double coe_1;
  for (i=0;i<=N-2;i++)
    {
      PiecewisePoly s;
      coe_0=y[i]-x[i]*((y[i+1]-y[i])/(x[i+1]-x[i]));
      s.coe.push_back(coe_0);
      coe_1=(y[i+1]-y[i])/(x[i+1]-x[i]);
      s.coe.push_back(coe_1);
      s.exp.push_back(0);
      s.exp.push_back(1);
      s.begin = x[i];
      s.end = x[i+1];
      answer.push_back(s);
    }
  return answer;
}

vector<PiecewisePoly> ppForm::S_23(int _condition,double _boundA,double _boundB)
{
 
  condition = _condition;
  boundA = _boundA;
  boundB = _boundB;
  vector<double> lambda;
  vector<double> mu;
  double lam;
  double miu;
  vector<vector<double>> table_dd;

  if(condition >3 || condition <=0)
    {
     cout << "Error Input!!" << endl;

     exit(0);

    }

  //计算lambda_i和mu_i并储存成向量
  
  for (int i=1;i<=N-2;i++) //lambda[0]=lambda2
    {
      lam=(x[i+1]-x[i])/(x[i+1]-x[i-1]);
      lambda.push_back(lam);
      miu=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
      mu.push_back(miu);
      }
  

  //计算插商表
  
  for (int i=0; i<N; i++)
    {
      vector<double> line;
      for (int j=0;j<N;j++)
	{
	  line.push_back(0);
	}
      table_dd.push_back(line);
    }
  
      for (int i=0;i<N;i++)
	{
	  if (i==0)
	    {
	      for (int j=0;j<N;j++)
		{
		  if (j==0)
		      table_dd[i][j]=y[i];
		  else 
		      table_dd[i][j]=0;
	        }
	    }
	  else
	    {
	      for (int j=0;j<N;j++)
		{
		  if (j==0)
		      table_dd[i][j]=y[i];
		  else if (j>i)
		      table_dd[i][j]=0;
		  else
		      table_dd[i][j]=(table_dd[i][j-1]-table_dd[i-1][j-1])/(x[i]-x[i-j]);
	        }
	     }
	}
  

  
  //计算D1样条
  if (condition == 1)
    {
      vector<double> bi;
      double bi2;
      bi2=3*mu[0]*table_dd[2][1]+3*lambda[0]*table_dd[1][1]-lambda[0]*boundA;
      bi.push_back(bi2);
      for (int i=1;i<=N-4;i++)
	{
	  double bii;
	  bii=3*mu[i]*table_dd[i+2][1]+3*lambda[i]*table_dd[i+1][1];
	  bi.push_back(bii);
	}
      double bil;
      bil=3*mu[N-3]*table_dd[N-1][1]+3*lambda[N-3]*table_dd[N-2][1]-mu[N-3]*boundB;
      bi.push_back(bil);
    

  //解三对角线性方程组
	vector<double>alpha(N-2);
	vector<double>beta(N-2);
	vector<double>Y(N-2);
	vector<double>m(N);
	m[0]=boundA;
	m[N-1]=boundB;
	alpha[0]=2;
	beta[0]=mu[0]/2;
	Y[0]=bi[0]/2;
	for (int i=1;i<N-2;i++)
	  {
	    alpha[i]=2-lambda[i]*beta[i-1];
	    if (i<N-3)
	      {
		beta[i]=mu[i]/alpha[i];
	      }
	    Y[i]=(bi[i]-lambda[i]*Y[i-1])/alpha[i];
	  }
	m[N-2]=Y[N-3];
        for (int i = N- 4; i >= 0; i--)
	  {
	    m[i+1] = Y[i] - beta[i] * m[i + 2];
          }
	
  //求出分段多项式
  
  for (int i=0;i<=N-2;i++)
    {
      PiecewisePoly s;
      double A;
      double B;
      double C;
      double coe_0;
      double coe_1;
      double coe_2;
      double coe_3;
      
      A=m[i];
      B=(table_dd[i+1][1]-m[i])/(x[i+1]-x[i]);
      C=(m[i+1]+m[i]-2*table_dd[i+1][1])/(x[i+1]*x[i+1]+x[i]*x[i]-2*x[i+1]*x[i]);
      
      coe_0=y[i]-x[i]*A+x[i]*x[i]*B-x[i]*x[i]*x[i+1]*C;
       s.coe.push_back(coe_0);
      coe_1=A-2*x[i]*B+(x[i]*x[i]+2*x[i]*x[i+1])*C;
      s.coe.push_back(coe_1);
      coe_2=B+(-x[i+1]-2*x[i])*C;
      s.coe.push_back(coe_2);
      coe_3=C;
      s.coe.push_back(coe_3);
      s.exp.push_back(0);
      s.exp.push_back(1);
      s.exp.push_back(2);
      s.exp.push_back(3);
      s.begin = x[i];
      s.end = x[i+1];
      answer.push_back(s);
    }
   return answer;

    }else if (condition == 2 || condition == 3)//计算D2样条
    {
      if (condition == 3)
	{
	  if (boundA != 0 || boundB != 0)
	    {
	      cout << "Error Input!!" << endl;
	      exit(0);
	    }
	}
      
      vector<double> bi;
      double bi2;
      bi2=6*table_dd[2][2]-mu[0]*boundA;
      bi.push_back(bi2);
      for (int i=1;i<=N-4;i++)
	{
	  double bii;
	  bii=6*table_dd[i+2][2];
	  bi.push_back(bii);
	}
      double bil;
      bil=6*table_dd[N-1][2]-lambda[N-3]*boundB;
      bi.push_back(bil);
      
    

  //解三对角线性方程组
      
	vector<double>alpha(N-2);
	vector<double>beta(N-2);
	vector<double>Y(N-2);
	vector<double>M(N);
	M[0]=boundA;
	M[N-1]=boundB;
	alpha[0]=2;
	beta[0]=lambda[0]/2;
	Y[0]=bi[0]/2;
	for (int i=1;i<N-2;i++)
	  {
	    alpha[i]=2-mu[i]*beta[i-1];
	    if (i<N-3)
	      {
		beta[i]=lambda[i]/alpha[i];
	      }
	    Y[i]=(bi[i]-mu[i]*Y[i-1])/alpha[i];
	  }
	M[N-2]=Y[N-3];
        for (int i = N- 4; i >= 0; i--)
	  {
	    M[i+1] = Y[i] - beta[i] * M[i + 2];
          }
      

  //求出分段多项式
	
  
  for (int i=0;i<=N-2;i++)
    {
      PiecewisePoly s;
      double A;
      double B;
      double C;
      double coe_0;
      double coe_1;
      double coe_2;
      double coe_3;
      
      
      A=table_dd[i+1][1]-((M[i+1]+2*M[i])*(x[i+1]-x[i]))/6;
      B=M[i]/2;
      C=((M[i+1]-M[i])/(x[i+1]-x[i]))/6;
      
      coe_0=y[i]-x[i]*A+x[i]*x[i]*B-x[i]*x[i]*x[i]*C;
      s.coe.push_back(coe_0);
       coe_1=A-2*x[i]*B+3*x[i]*x[i]*C;
      s.coe.push_back(coe_1);
      coe_2=B-3*x[i]*C;
      s.coe.push_back(coe_2);
      coe_3=C;
      s.coe.push_back(coe_3);
      s.exp.push_back(0);
      s.exp.push_back(1);
      s.exp.push_back(2);
      s.exp.push_back(3);
      s.begin = x[i];
      s.end = x[i+1];
      answer.push_back(s);
    }
      
   return answer;
    }
  return answer;

  }

class Bspline : public Solve_Base
{
 public:
  Bspline(int _N)
    :Solve_Base(_N)
  {
    
  }
  Bspline(int _N,vector<double> _x,vector<double> _y)
    :Solve_Base(_N,_x,_y)
  {
    N=_N;
    for (int i=0;i<=N-1;i++)
      {
	x.push_back(_x[i]);
	y.push_back(_y[i]);
      }
  }

  
  vector<PiecewisePoly> S_01();
  vector<PiecewisePoly> S_23(int,double,double);
  vector<PiecewisePoly> S_12(double,double);
  int condition;
  double boundA;
  double boundB;
  
};

vector<PiecewisePoly> Bspline::S_01()
{
  int i;
  double coe_0;
  double coe_1;
  for (i=0;i<=N-2;i++)
    {
      PiecewisePoly s;
      coe_0=(y[i]*x[i+1]-y[i+1]*x[i])/(x[i+1]-x[i]);
      s.coe.push_back(coe_0);
      coe_1=(-y[i]+y[i+1])/(x[i+1]-x[i]);
      s.coe.push_back(coe_1);
      s.exp.push_back(0);
      s.exp.push_back(1);
      s.begin = x[i];
      s.end = x[i+1];
      answer.push_back(s);
    }
  return answer;
}

vector<PiecewisePoly> Bspline::S_23(int _condition,double _boundA,double _boundB)
{
  condition = _condition;
  boundA = _boundA;
  boundB = _boundB;
  
  if(condition >3 || condition <=0)
    {
     cout << "Error Input!!" << endl;

     exit(0);

    }
  
  //计算D1样条
  if (condition == 1)
    {
      vector<double> bi;
      double bi1;
      bi1=3*y[0]+boundA;
      bi.push_back(bi1);
      for (int i=1;i<=N-2;i++)
	{
	  double bii;
	  bii=6*y[i];
	  bi.push_back(bii);
	}
      double bil;
      bil=3*y[N-1]-boundB;
      bi.push_back(bil);
    

  //解三对角线性方程组
	vector<double>alpha(N);
	vector<double>beta(N);
	vector<double>Y(N);
	vector<double>a(N+2);
	alpha[0]=2;
	beta[0]=0.5;
	Y[0]=bi[0]/2;
	for (int i=1;i<N;i++)
	  {
	    if (i==N-1)
	      {
		alpha[i]=2-beta[i-1];
	      }
	    else
	      {
		alpha[i]=4-beta[i-1];
	      }
	    
	    if (i<N-1)
	      {
		beta[i]=1/alpha[i];
	      }
	    Y[i]=(bi[i]-Y[i-1])/alpha[i];
	  }
	a[N]=Y[N-1];
        for (int i = N-2; i >= 0; i--)
	  {
	    a[i+1] = Y[i] - beta[i] * a[i + 2];
          }
	a[0]=a[2]-2*boundA;
	a[N+1]=a[N-1]+2*boundB;

  //求出分段多项式
  
  for (int i=1;i<=N-1;i++)
    {
      PiecewisePoly s;
      double coe_0;
      double coe_1;
      double coe_2;
      double coe_3;
      double i_=i;
      
      coe_0=(a[i-1]/6)*pow(i_+1,3)+((2*a[i])/3)-0.5*a[i]*(i_*i_*i_+2*i_*i_)+((2*a[i+1])/3)+0.5*a[i+1]*(i_-1)*(i_+1)*(i_+1)-(a[i+2]/6)*pow(i_,3);
      s.coe.push_back(coe_0);
      coe_1=(-a[i-1]/2)*(i_+1)*(i_+1)+(a[i]/2)*(3*i_*i_+4*i_)-(a[i+1]/2)*(3*i_*i_+2*i_-1)+(a[i+2]/2)*(i_)*(i_);
      s.coe.push_back(coe_1);
      coe_2=(a[i-1]/2)*(i_+1)-(a[i]/2)*(3*i_+2)+(a[i+1]/2)*(3*i_+1)-(a[i+2]/2)*(i_);
      s.coe.push_back(coe_2);
      coe_3=-(a[i-1]/6)+a[i]/2-a[i+1]/2+a[i+2]/6;
      s.coe.push_back(coe_3);
      s.exp.push_back(0);
      s.exp.push_back(1);
      s.exp.push_back(2);
      s.exp.push_back(3);
      s.begin = x[i-1];
      s.end = x[i];
      answer.push_back(s);
    }
   return answer;

    }else if (condition == 2 || condition == 3)//计算D2样条
    {
      if (condition == 3)
	{
	  if (boundA != 0 || boundB != 0)
	    {
	      cout << "Error Input!!" << endl;
	      exit(0);
	    }
	}
      double wei;
      wei=x[0]-1;
      vector<double> bi;
      double bi1;
      bi1=6*y[0]-boundA;
      bi.push_back(bi1);
      for (int i=1;i<=N-2;i++)
	{
	  double bii;
	  bii=6*y[i];
	  bi.push_back(bii);
	}
      double bil;
      bil=6*y[N-1]-boundB;
      bi.push_back(bil);
    

  //解三对角线性方程组

	vector<double>alpha(N);
	vector<double>beta(N);
	vector<double>Y(N);
	vector<double>a(N+2);
	alpha[0]=6;
	beta[0]=0;
	Y[0]=bi[0]/6;
	for (int i=1;i<N;i++)
	  {
	    if (i==N-1)
	      {
		alpha[i]=6;
	      }
	    else
	      {
		alpha[i]=4-beta[i-1];
	      }
	    
	    if (i<N-1)
	      {
		beta[i]=1/alpha[i];
	      }
	    if (i==N-1)
	      {
		Y[i]=(bi[i])/alpha[i];
	      }
	    else
	      {
		Y[i]=(bi[i]-Y[i-1])/alpha[i];
	      }
	    
	  }
	a[N]=Y[N-1];
        for (int i = N-2; i >= 0; i--)
	  {
	    a[i+1] = Y[i] - beta[i] * a[i + 2];
          }
	a[0]=boundA+2*a[1]-a[2];
	a[N+1]=boundB+2*a[N]-a[N-1];

	
  //求出分段多项式
  
  for (int i=1;i<=N-1;i++)
    {

      PiecewisePoly s;
      double coe_0;
      double coe_1;
      double coe_2;
      double coe_3;
      
      coe_0=(a[i-1]/6)*pow(i+1,3)+((2*a[i])/3)-0.5*a[i]*(i*i*i+2*i*i)+((2*a[i+1])/3)+0.5*a[i+1]*(i-1)*(i+1)*(i+1)-(a[i+2]/6)*pow(i,3);
      s.coe.push_back(coe_0);
      coe_1=(-a[i-1]/2)*(i+1)*(i+1)+(a[i]/2)*(3*i*i+4*i)-(a[i+1]/2)*(3*i*i+2*i-1)+(a[i+2]/2)*(i)*(i);
      s.coe.push_back(coe_1);
      coe_2=(a[i-1]/2)*(i+1)-(a[i]/2)*(3*i+2)+(a[i+1]/2)*(3*i+1)-(a[i+2]/2)*(i);
      s.coe.push_back(coe_2);
      coe_3=-(a[i-1]/6)+a[i]/2-a[i+1]/2+a[i+2]/6;
      s.coe.push_back(coe_3);
      s.exp.push_back(0);
      s.exp.push_back(1);
      s.exp.push_back(2);
      s.exp.push_back(3);
      s.begin = x[i-1];
      s.end = x[i];
      answer.push_back(s);
      
    }
   return answer;
    }
  return answer;

  }

vector<PiecewisePoly> Bspline::S_12(double _boundA,double _boundB)
{
  boundA = _boundA;
  boundB = _boundB;


      vector<double> bi;
      double bi1;
      bi1=8*y[0]-2*boundA;
      bi.push_back(bi1);
      for (int i=1;i<=N-2;i++)
	{
	  double bii;
	  bii=8*y[i];
	  bi.push_back(bii);
	}
      double bil;
      bil=8*y[N-1]-2*boundB;
      bi.push_back(bil);
    

  //解三对角线性方程组
	vector<double>alpha(N);
	vector<double>beta(N);
	vector<double>Y(N);
	vector<double>a(N+2);
	alpha[0]=5;
	beta[0]=0.2;
	Y[0]=bi[0]/5;
	for (int i=1;i<N;i++)
	  {
	    if (i==N-1)
	      {
		alpha[i]=5-beta[i-1];
	      }
	    else
	      {
		alpha[i]=6-beta[i-1];
	      }
	    
	    if (i<N-1)
	      {
		beta[i]=1/alpha[i];
	      }
	    Y[i]=(bi[i]-Y[i-1])/alpha[i];
	  }
	a[N]=Y[N-1];
        for (int i = N-2; i >= 0; i--)
	  {
	    a[i+1] = Y[i] - beta[i] * a[i + 2];
          }
	a[0]=2*boundA-a[1];
	a[N+1]=-a[N]+2*boundB;

  //求出分段多项式
  
  for (int i=1;i<=N;i++)
    {
      PiecewisePoly s;
      double coe_0;
      double coe_1;
      double coe_2;
      double i_=i;
      
      coe_0=(a[i-1]/2)*pow(i_+1,2)+0.75*a[i]-a[i]*(i_*i_+0.25+i_)+0.5*a[i+1]*(i_)*(i_);
      s.coe.push_back(coe_0);
      coe_1=-a[i-1]*(i_+1)+a[i]*(2*i_+1)-a[i+1]*(i_);
      s.coe.push_back(coe_1);
      coe_2=(a[i-1]/2)-a[i]+(a[i+1]/2);
      s.coe.push_back(coe_2);
      s.exp.push_back(0);
      s.exp.push_back(1);
      s.exp.push_back(2);
      s.begin = x[i-1];
      if (i!=N)
	{
	  s.end = x[i];
	}
      else
	{
	  s.end = s.begin+1;
	}
      
      answer.push_back(s);
    }
   return answer;
    }



class Curve : public Solve_Base
{
 public:
  Curve(int _N,vector<double> _x,vector<double> _y )
    :Solve_Base(_N,_x,_y)
  {
    vector<double> t_(N);
    for (int i=0;i<N;i++)
      {
	if (i==0)
	  {
	    t_[i]=0;
	  }
	else
	  {
	    t_[i]=t_[i-1]+pow(pow(x[i]-x[i-1],2)+pow(y[i]-y[i-1],2),1/2);
	  }
      }
    
    vector<double> x_real(N+1);
    for (int i=0;i<N;i++)
      {
	x_real[i]=x[i];
      }
    x_real[N]=x[0];

    vector<double> y_real(N+1);
    for (int i=0;i<N;i++)
      {
	y_real[i]=y[i];
      }
    y_real[N]=y[0];
    for (int i=0;i<N;i++)
      {
	t.push_back(t_[i]);
      }
    double tt=t_[N-1]+pow(pow(x_real[N]-x[N-1],2)+pow(y[N]-y[N-1],2),1/2);
    t.push_back(tt);
    
    ppForm cur(N+1,t,x_real);
    answer= cur.S_23(3,0,0);

    ppForm ty(N+1,t,y_real);
    answer_=ty.S_23(3,0,0);
    
  }

};
